Post-Newtonian Cosmological Dynamics 
Lagrangian coordinates 



Sabino Matarrese & David Terranova 

Dipartimento di Fisica "Galileo Galilei", 
Universita di Padova, via Marzolo 8, 35131 Padova, Italy 

February 5, 2008 

Abstract 

We study the non-linear dynamics of self-gravitating irrotational dust in a general 
relativistic framework, using synchronous and comoving (i.e. Lagrangian) coordinates. 
All the equations are written in terms of a single tensor variable, the metric tensor of 
the spatial sections orthogonal to the fluid flow. This treatment allows an unambiguous 
expansion in inverse (even) powers of the speed of light. To lowest order, the Newtonian 
approximation - in Lagrangian form - is derived and written in a transparent way; the 
corresponding Lagrangian Newtonian metric is obtained. Post-Newtonian corrections 
are then derived and their physical meaning clarified. A number of results are obtained: 
i) the master equation of Lagrangian Newtonian dynamics, the Raychaudhuri equation, 
can be interpreted as an equation for the evolution of the Lagrangian-to-Eulerian Jaco- 
bian matrix, complemented by the irrotationality constraint; ii) the Lagrangian spatial 
metric reduces, in the Newtonian limit, to that of Euclidean 3-space written in time- 
dependent curvilinear coordinates, with non-vanishing Christoffel symbols, but vanishing 
spatial curvature (a particular example of it is given within the Zel'dovich approxima- 
tion); in) a Lagrangian version of the Bernoulli equation for the evolution of the "velocity 
potential" is obtained, iv) The Newtonian and post-Newtonian content of the electric 
and magnetic parts of the Weyl tensor is clarified, v) At the Post-Newtonian level, an 
exact and general formula is derived for gravitational-wave emission from non-linear 
cosmological perturbations; vi) a straightforward application to the anisotropic collapse 
of homogeneous ellipsoids shows that the ratio of these post-Newtonian terms to the 
Newtonian ones tends to diverge at least like the mass density, vii) It is argued that a 
stochastic gravitational-wave background is produced by non-linear cosmic structures, 
with present-day closure density Q gw ~ 10~ 5 - 10~ 6 on 1 - 10 Mpc scales. 

Key words: gravitation - hydrodynamics - instabilities - cosmology: theory - large- 
scale structure of Universe. 



1 Introduction 



The gravitational instability of collisionless matter in a cosmological framework is usually stud- 
ied within the Newtonian approximation, which basically consists in neglecting terms whose 
order is higher than the first in metric perturbations around a matter-dominated Friedmann- 
Robertson- Walker (FRW) background, while keeping non-linear density perturbations. This 
approximation is usually thought to produce accurate results in a wide spectrum of cosmolog- 
ical scales, namely on scales much larger than the Schwarzschild radius of collapsing bodies 
and much smaller than the Hubble horizon scale, where the peculiar gravitational potential 
ip g , divided by the square of the speed of light c 2 to obtain a dimensionless quantity, keeps 
much less than unity, while the peculiar matter flow never becomes relativistic. To be more 
specific, the Newtonian approximation consists in perturbing only the time-time component 
of the FRW metric tensor by an amount 2ip g /c 2 , where ip g is related to the matter density 
fluctuation 5 via the cosmological Poisson equation, 

V 2 ^(x, t) = 47rGa 2 (t) Qb (t)5(x, t) , (1) 

where g b is the background matter density and a(t) the appropriate FRW scale-factor; the 
Laplacian operator V 2 has been used here with its standard meaning of Euclidean space. 
The fluid dynamics is then usually studied in Eulerian coordinates by accounting for mass 
conservation and using the cosmological version of the Euler equation for a self-gravitating 
pressureless fluid - as long as the flow is in the laminar regime - to close the system. To 
motivate the use of this "hybrid approximation", which deals with perturbations of the matter 
and the geometry at a different perturbative order, one can either formally expand the correct 
equations of General Relativity (GR) in inverse powers of the speed of light (e.g. Weinberg 
1972), or simply notice that the peculiar gravitational potential is strongly suppressed with 
respect to the matter perturbation by the square of the ratio of the perturbation scale A to 
the Hubble radius r# = cH" 1 (H being the Hubble constant): y? 3 /c 2 ~ 5 (X/th) 2 - 

Such a simplified approach, however, already fails in producing an accurate description 
of the trajectories of relativistic particles, such as photons. Neglecting the relativistic per- 
turbation of the space-space components of the metric, which in the so-called longitudinal 
gauge is just —2(p g /c 2 , would imply a mistake by a factor of two in well-known effects such 
as the Rees-Sciama (1968) and gravitational lensing (e.g. Schneider, Ehlers & Falco 1992), as 
it would be easy to see, by looking at the solution of the eikonal equation. In other words, 
the level of accuracy not only depends on the peculiar velocity of the matter producing the 
spacetime curvature, but also on the nature of the particles carrying the signal to the observer. 
Said this way, it may appear that the only relativistic correction required to the usual Eulerian 
Newtonian picture is that of writing the metric tensor in the revised, "weak field", form (e.g. 
Peebles 1993) 

ds 2 = - (l + ^) c 2 dt 2 + a 2 (t) (l - ^) dl 2 . (2) 

However, as we are going to show, this is not the whole story. It is well-known in fact that 
the gravitational instability of aspherical perturbations (which is the generic case) leads to the 
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formation of very anisotropic structures whenever pressure gradients can be neglected (Lynden- 
Bell 1962; Lin, Mestel & Shu 1965; Zel'dovich 1970; Icke 1973; White & Silk 1979; Shandarin 
et al. 1995). Matter first flows in almost two-dimensional structures called pancakes, which 
then merge and fragment to eventually form one-dimensional filaments and point-like clumps. 
During the process of pancake formation the matter density, the shear and the tidal field 
formally become infinite along evanescent two-dimensional configurations corresponding to 
caustics; after this event a number of highly non-linear phenomena, such as vorticity generation 
by multi-streaming, merging, tidal disruption and fragmentation, occur. Most of the patology 
of the caustic formation process, such as the local divergence of the density, shear and tide, and 
the formation of multi-stream regions, are just an artifact of extrapolating the pressureless fluid 
approximation beyond the point at which pressure gradients and viscosity become important. 
In spite of these limitations, however, it is generally believed that the general anisotropy of 
the collapse configurations, either pancakes or filaments, is a generic feature of cosmological 
structures originated through gravitational instability, which would survive even in the presence 
of a collisional component. 

This simple observation already shows the inadequacy of the standard Newtonian paradigm. 
According to it, the lowest scale at which the approximation can be reasonably applied is set 
by the amplitude of the gravitational potential and is given by the Schwarzschild radius of 
the collapsing body, which is negligibly small for any relevant cosmological mass scale. What 
is completely missing in this criterion is the role of the shear which causes the presence of 
non-scalar contributions to the metric perturbations. A non-vanishing shear component is in 
fact an unavoidable feature of realistic cosmological perturbations and affects the dynamics in 
(at least) three ways, all related to non-local effects, i.e. to the interaction of a given fluid 
element with the environment. 

First, at the lowest perturbative order the shear is related to the tidal field generated by 
the surrounding material by a simple proportionality law (because of this linear coincidence, 
in much of the literature "shear" and "tide" are used as synonima). This sort of non-locality, 
however, is coded in the initial conditions of each fluid-element through a Coulomb-like in- 
teraction with arbitrarily distant matter. Because of its link with the initial data of each fluid 
element one can consider it as a local property. The later modification of these shear and tidal 
fields is one of the consequences of the non-linear evolution. 

Second, it is related to a dynamical tidal induction: the modification of the environment 
forces the fluid element to modify its shape and density. In Newtonian gravity, this is an 
action- at- a- distance effect, which starts to manifest itself in second-order perturbation theory 
as an inverse-Laplacian contribution to the velocity potential (e.g. Catelan et al. 1995, and 
references therein). 

Third, and most important here, a non-vanishing shear field leads to the generation of 
a traceless and divergenceless metric perturbation which can be understood as gravitational 
radiation emitted by non-linear perturbations. This contribution to the metric perturbations 
is statistically small on cosmologically interesting scales, but it becomes relevant whenever 
anisotropic (with the only exception of exactly one-dimensional) collapse takes place. In the 
Lagrangian picture considered here, such an effect already arises at the post-Newtonian (PN) 
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level. 

Note that the two latter effects are only detected if one allows for non-scalar perturbations 
in physical quantities. Contrary to a widespread belief, in fact, the choice of scalar pertur- 
bations in the initial conditions is not enough to prevent tensor modes to arise beyond the 
linear regime in a GR treatment. Truly tensor perturbations are dynamically generated by the 
gravitational instability of initially scalar perturbations, independently of the initial presence 
of gravitational waves. 

This point is very clearly displayed in the GR Lagrangian second-order perturbative ap- 
proach. The pioneering work in this field is by Tomita, who, back in 1967, calculated the grav- 
itational waves emitted by non-linearly evolving scalar perturbations in an Einstein-de Sitter 
background, in the synchronous gauge (Tomita 1967). Matarrese, Pantano & Saez (1994a,b) 
obtained an equivalent result but with a different formalism in comoving and synchronous 
coordinates. According to these calculations, a traceless and divergenceless contribution to 
the spatial metric in the synchronous gauge, it 01 * [greek indices label Lagrangian spatial coordi- 
nates, while capital latin letters will label Eulerian space; lower-case latin indices will be used 
for spacetime coordinates], is produced, which, with growing mode initial conditions, obeys 
the inhomogeneous wave-equation 

+ ^ ~ V «N = -42 v « s " ' (3) 

where r oc t 1 ^ 3 is the conformal time. The non-linear source tensor S a p is given in terms of 
the initial peculiar gravitational potential, <po = <p g (to), by 

S a p = S" p Vfa + VC/j + 2 (<Po % Vfco - VXf^) , (4) 

with 

V^o = -^((V^o) 2 -^o,>o;,) , (5) 

where spatial gradients, indicated by greek indices after a comma, are with respect to the 
Lagrangian coordinates q a and indices are raised by the Kronecker symbol; finally V 2 , is just 
the standard (i.e. Euclidean) Laplacian in Lagrangian coordinates. To get from the above 
formula a form which can be compared with the standard Newtonian interpretation, we can 
expand n a p in powers of 1/c 2 (as we will see below, the absence of odd powers of the speed of 
light is a characteristic feature of the Lagrangian approach), ir^ = 7T^ a ^ + ^^ PN ^ a p + C(^)- 
To zeroth order one obtains the Newtonian term ir^^g = Q, which includes a non-local 

and non-causal contribution to the shear tensor through derivatives of the potential ipo defined 
above. The meaning of this contribution has been discussed by Matarrese, Pantano & Saez 
(1994a,b), who obtained it by looking at perturbation scales much smaller than the Hubble 
radius; it represents the "relic" of a causal signal which, on sub-horizon scales, appears as 
an instantaneous Newtonian feature. To first order in 1/c 2 one then gets a PN contribution, 

1 Actually, this expression is determined up to a harmonic divergenceless and traceless tensor, which can be 
set to zero if we require consistency with the standard Newtonian second-order results. 
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Vg7r( pjv ) a = ^-S 01 ^, once again the causality of this gravitational-wave signal is lost because 
of the 1/c 2 expansion. In the formalism of (Matarrese, Pantano & Saez 1994a,b) this PN term 
would be detected as a sub-leading contribution for perturbation scales much smaller than the 
Hubble radius. The close relation between the two approximation schemes - inverse powers of 
the speed of light and powers of the ratio of the perturbation to the horizon scale - also helps 
in better understanding the actual physical meaning of the 1/c 2 expansion in the Lagrangian 
picture. 

The latter PN effect will be recovered in Section 4 without any restriction to a second-order 
perturbation treatment. A heuristic estimate of the amplitude of this effect in the frame of 
current scenarios of cosmological structure formation is reported in Section 5. One can also 
speculate on the possibility to detect the resulting stochastic gravitational-wave background, 
e.g., through the secondary anisotropy it would induce on the Cosmic Microwave Background 
(CMB). 

There is an intriguing aspect of the above expression for the tensor modes in the PN limit, 
namely V 2 7r( PAr) a = ^-S%. Our explicit PN result of Section 4, if further approximated to a 
perturbative second-order expansion around the FRW background, gives V^K^ PN ^ a ^ = ^§S a pi 
i.e. a factor of 6 smaller! Should one conclude that the 1/c 2 expansion and that in the 
amplitude of the perturbations around FRW do not commute? The explanation is actually 
much simpler: the splitting of the various perturbation modes into scalars, vectors ans tensors is 
obviously background-dependent; in particular, our PN tensor modes are defined with respect 
to a non-perturbative Newtonian background, while the tensor modes obtained in second- 
order perturbation theory are defined with respect to the usual FRW solution. When the 
PN expressions for the various geometric modes are further expanded to second order in 
perturbation theory and then compared to those obtained through a second-order followed by 
a 1/c 2 expansion the results do not generally coincide (because the Newtonian background itself 
contains perturbations of the FRW metric), but their sum, i.e. the overall metric perturbation, 
does. This completely accounts for the different numerical factor. 

Finally, Eq.(3) allows to understand another important point: the complete insensitivity of 
the Newtonian approximation to the possible presence of free gravitational waves in the initial 
conditions, such as those produced by quantum effects in the early universe. These initial 
tensor modes, corresponding to solutions of the homogeneous equation associated to Eq.(3), 
would reduce to harmonic, transverse and traceless metric perturbations in the Newtonian 
limit, having no effect on physical quantities (they are gauge modes from the point of view of 
the Newtonian equations). 

The reader at this point may be confused by the continuous interchange of Newtonian and 
PN concepts. However, this will appear unavoidable once one realizes that, as in any pertur- 
bative treatment (the perturbation parameter here being formally 1/c 2 ), there are equations 
which mix different perturbation orders. So, the PN equations will have Newtonian sources, 
or read the other way around, there are Newtonian effects which are produced by PN sources. 
This point has been definitely clarified in a fundamental paper by Kofman & Pogosyan (1995), 
who showed how the Newtonian "electric" tidal field E a p evolves in time according to a PN 
equation, so that the circulation of the PN "magnetic" Weyl tensor H 01 *, happens to be respon- 
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sible for the Newtonian non-local "tidal induction". Bertschinger & Hamilton (1994) gave a 
different interpretation of the same effect. 

Recently a number of different approaches to relativistic effects in the non-linear dynamics 
of cosmological perturbations have been proposed. Matarrese, Pantano & Saez (1993) pro- 
posed an algorithm based on neglecting the magnetic part of the Weyl tensor in the dynamics, 
obtaining strictly local fluid-flow evolution equations, i.e. the so-called "silent universe". Us- 
ing this formalism Bruni, Matarrese & Pantano (1995a) studied the asymptotic behaviour of 
the system, both for collapse and expansion, showing, in particular, that this kind of local 
dynamics generically leads to spindle singularities for collapsing fluid elements, thereby con- 
firming the results of a previous analysis by Bertschinger & Jain (1994). This formalism, 
however, cannot be applied to cosmological structure formation inside the horizon, where the 
non-local tidal induction cannot be neglected, i.e. the magnetic Weyl tensor H a g is non-zero, 
with the exception of highly specific initial configurations (Matarrese, Pantano & Saez 1994a; 
Bertschinger & Jain 1994). Rather, it is probably related to the non-linear dynamics of an ir- 
rotational fluid outside the (local) horizon (Matarrese, Pantano & Saez 1994a,b). One possible 
application (Bruni, Matarrese & Pantano 1995b), is in fact connected to the so-called Cosmic 
No-hair Theorem (e.g. Hawking & Moss 1982), i.e. to the conjecture that expanding patches of 
an initially inhomogeneous and anisotropic universe asymptotically tend to almost FRW solu- 
tions, thanks to the action of a cosmological constant-like term. The self-consistency of these 
"silent universe" models has been recently demonstrated by Lesame, Dunsby & Ellis (1995), 
extending an earlier analysis by Barnes & Rowlingson (1989). Lesame, Ellis & Dunsby (1996) 
discussed the role of the divergence of the magnetic Weyl tensoif], whose presence reflects the 
fact that the shear and the electric tide generally have a different eigenframe. A local-tide ap- 
proximation for the non-linear evolution of collisionless matter, which tries to overcome some 
limitations of the Zel'dovich approximation (Zel'dovich 1970), has been recently proposed by 
Hui & Bertschinger (1995). 

In this work we will follow the more "conservative" approach of expanding the Einstein and 
continuity equations in inverse powers of the speed of light, which will then define a Newtonian 
limit and, at the next order, post-Newtonian corrections. The newer aspect of our approach 
is the choice of gauge: we use synchronous and comoving coordinates, because of which our 
approach can be legitimately called a Lagrangian one. Thanks to this choice, the dynamical 
variables involved are quite different to the standard ones; the gravitational potential, for 
instance, never appears explicitly in our expansion. 

Various approaches have been proposed in the literature, which are somehow related to 
the present one. A PN approximation has been followed by Futamase (1988, 1989, 1991) to 
describe the dynamics of a clumpy universe; he however used non-comoving coordinates and 
focused his analysis on applications related to the so-called averaging problem in cosmology 
(e.g. Ellis 1984). Tomita (1988, 1991) also used non-comoving coordinates in a PN approach 
to cosmological perturbations. Shibata & Asada (1995) recently developed a PN approach to 
cosmological perturbations, but they also used non-comoving coordinates. Kasai (1995) [see 

2 Actually, a non-zero divH only appears as a third-order effect in the amplitude of perturbations around 
FRW, unless the initial conditions contain a mixture of growing and decaying modes. 
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also (Kasai 1992, 1993)] analyzed the non-linear dynamics of dust in the synchronous and 
comoving gauge; his approximation methods are however largely different. Finally, in a series 
of papers, based on the Hamilton- Jacobi approach (Croudace et al. 1994; Parry et al. 1994; 
Salopek, Stewart & Croudace 1994) a new approximation technique has been developed, which 
relies on an expansion in higher and higher gradients of an initial perturbation "seed" . In spite 
of its elegance and generality, however, this approximation scheme is by construction unable 
to reproduce the non-local aspects of the gravitational instability on sub-horizon scales; more 
specifically, terms containing the inverse of the Laplacian operator, which are unavoidable in 
the Newtonian limit, would formally require an infinite series of terms in a gradient expansion. 

To help the reader from not being too much confused by the various perturbative techniques 
adopted in this work, we anticipate that our calculations contain, in different parts, three 
different kinds of expansion. 

First, the entire paper is mostly based on an expansion in inverse even powers of the speed 
of light: the lowest order - or background - solution, in this case, describes the so-called 
Newtonian approximation in an expanding universe. Although in general we do not know 
the explicit form of the Newtonian background, we can safely assume it exists and use it to 
derive the next order terms. The result of first-order perturbation theory is then called post- 
Newtonian (PN); the second order would be the post-post-Newtonian (PPN) approximation. 
The range of application of this perturbative method has been already discussed above. Going 
to higher and higher orders would generally lead to a more accurate description of the system, 
account for some new relativistic effects, such as the generation of gravitational waves, and 
possibly allow for an extension of the range of scales to which the formalism can be applied. 

Second, we will also use the most standard cosmological perturbation theory (see, e.g., 
Kodama & Sasaki 1984, and references therein; for a pressureless medium, see also Hwang 
1994), which is basically an expansion in powers of the amplitude of the perturbations around 
a background, homogeneous and isotropic, FRW solution. The first-order, or "linear", terms 
of the expansion are given in Section 2.4. No second-order calculations, in this sense, will 
be presented here, with the exception of Eqs.(3) - (5) reported above. The works by Tomita 
(1967), Matarrese, Pantano & Saez (1994a,b) and - in some respect - Pyne & Carroll (1995) 
follow precisely this perturbative approach up to second order and within GR. The range of 
applicability of this second perturbation technique is that of small fluctuations around a FRW 
background, but with no extra limitations on scale. Going to higher and higher orders here 
generally helps to follow the gravitational instability process on a longer time-scale and to 
account for new non-linear and non-local phenomena. 

Third, there is another meaning of "perturbation theory" in Lagrangian coordinates, which 
is frequently used in the cosmological literature (e.g. Buchert 1995, and references therein). 
This refers to an expansion, within Newtonian gravity, in powers of the displacement vector 
from Lagrangian to Eulerian coordinates (Buchert 1989; Moutarde et al. 1991; Bouchet et 
al. 1992; Buchert 1992; Catelan 1995), the background being once more represented by the 
FRW models. The linear result is the so-called Zel'dovich approximation (see Section 3.2 
below), while the second order terms are either called "second-order Lagrangian" or "post- 
Zel'dovich" (e.g. Munshi, Sahni & Starobinsky 1994). The peculiarity of this treatment, 
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at any order, is that, while the displacement vector is calculated from the equations at the 
required perturbative order, all the other dynamical variables, such as mass density, shear and 
so on, are calculated exactly from their non-perturbative definition. What comes out is a fully 
non-linear description of the system, which, though not being generally correct, "mimics" the 
true non-linear behaviour. This perturbation treatment basically exploits the advantages of 
the Lagrangian picture, leading, in particular, to a more accurate description of high density 
regions. Its limitations are generally set by the emerging of caustic singularities, besides those 
deriving from the underlying Newtonian approximation. A similar, Zel'dovich-like, approach 
can also be followed within GR [some progress in this direction has been recently made by 
Kasai (1995)]; this will be however the subject of a future investigation. 

The plan of the paper is as follows. In Section 2 we introduce the GR Lagrangian formalism. 
Although we do not use the whole machinery of the ADM approach (Arnowitt, Deser & Misner 
1962), some of the language is the same; in particular, the clear distinction between constraint 
and evolution equations plays a key role also in our work. Section 3 deals with the Newtonian 
limit of the GR equations in Lagrangian coordinates and gives a number of formal applications 
of the approach. Section 4 is instead devoted to the post-Newtonian limit of the GR equations. 
In particular, we discuss the dynamical role of gravitational waves generated by non-linear 
cosmic structures. It is shown that, during the collapse of a non-spherical (and non-planar) 
perturbation, PN tensor modes are produced, so as to give a dominant contribution near the 
collapse time. Conclusions are drawn in Section 5, which also contains a qualitative discussion 
on the amplitude of the PN gravitational-wave modes, as well as some speculations on their 
possible detectability. 

2 Relativistic dynamics of irrotational dust in the La- 
grangian picture 

In this section we will derive the equations governing the evolution of an irrotational fluid of 
dust (i.e. p = u = 0) in a synchronous and comoving system of coordinates (actually the 
possibility of making these two gauge choices simultaneously is a peculiarity of irrotational 
dust, which holds at any time, i.e. also beyond the linear regime). The starting point will be 
the Einstein equations R a b — \g a bR = ^rT a b, with R a t, the Ricci tensor, and the continuity 
equation T ab . a = for the matter stress-energy tensor T ah = gc 2 u a u b , where g is the mass 
density and u a the fluid four-velocity (normalized to u a u a = —1); a semicolon denotes covariant 
differentiation. The line-element reads 

ds 2 = -c 2 dt 2 + /i Q/3 (q, t)dq a dq P . (6) 

The fluid four-velocity in comoving coordinates is u a = (1, 0, 0, 0). A fundamental quantity of 
our analysis will be the velocity-gradient tensor, which is purely spatial, 

9^ = cu^ = IhTkyf, , (7) 
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where a dot denotes partial differentiation with respect to the proper time t. The tensor 0°^ 
represents the extrinsic curvature of the spatial hypersurfaces orthogonal to u a . 

Thanks to the spacetime splitting obtained in our frame, the 10 Einstein equations can be 
immediately divided into 4 constraints and 6 evolution equations. The time-time component 
of the Einstein equations is the so-called energy constraint of the ADM approach, which reads 

e 2 - 0^0^ + c 2 ®r = iqtxGq , (8) 

where the volume-expansion scalar is just the trace of the velocity-gradient tensor, 

is the trace of the three-dimensional Ricci curvature, ^R a p, of the spatial hypersurfaces of 

constant time. 

The space-time components give the momentum constraint, 

©V = ,/3 • (9) 

Finally, the space-space components represent the only truly evolution equations, i.e. those 
which contain second-order time derivatives of the metric tensor. They indeed govern the 
evolution of the extrinsic curvature tensor and read 

+ 00^ + c 2 wk* p = AnGQ5 a p . (10) 

Taking the trace of the last equation and combining it with the energy constraint, we obtain 
the Raychaudhuri equation (Raychaudhuri 1957), 

+ 0^0^ + 4^ = 0. (11) 

Mass conservation is provided by the equation 

q=-Gq. (12) 

Given that = \h a ^h ia = dihih 1 ^ 2 ) / 'dt, where h = det h a p, we can write the solution of 
this equation in the form 

^(q,t) = ^o(q)[^(q,t)Ao(q)]" 1/2 . (13) 

Here and in what follows quantities with a subscript are meant to be evaluated at some 
initial time to- 

Finally, let us introduce the so-called electric and magnetic parts of the Weyl tensor, which 
are both symmetric, flow-orthogonal and traceless. They read, respectively, 

E% = 1*°, (©^ - 2 ) + 00^ - Q^Q\ + c 2 ( - l -5% (%) (14) 

and 

H% = \h Plt fr l5 ® a r , S + V^Q^s) , (15) 

where rf^ = h~ l l 2 e al51 is the three-dimensional, completely anti-symmetric Levi-Civita ten- 
sor relative to the spatial metric h a p and e Q/37 is such that e 123 = 1, etc... . 
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Notice that, while the definition of the electric tide E a p is completely fixed, because of 
its well-known Newtonian limit, the magnetic tensor field has no straightforward Newtonian 
counterpart, and can be therefore defined up to arbitrary powers of the speed of light. The 
definition we are adopting here is the most straightforward one; it is such that no explicit 
powers of c appear in Eq.(15), which means that its physical dimensions are 1/c those of E a g. 
This choice can be motivated in analogy with electrodynamics, where the magnetic vector field 
is also scaled by 1/c with respect to the electric one. We will come back later, in Section 2.2 
and Section 4.1, to the consequences of this choice. 
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2.1 Conformal rescaling and FRW background subtraction 

With the purpose of studying gravitational instability in a FRW background, it is convenient 
to factor out the homogeneous and isotropic solutions of the above equations. To this aim we 
also perform a conformal rescaling of the metric with conformal factor a(t), the scale-factor of 
FRW models, and change the time variable to the conformal time r, defined by dr = dt/a{t). 
The line-element is then written in the form 

ds 2 = a 2 {r)[ - c 2 dr 2 + 7a/3 (q, r)dq a dq^} , (16) 

where a 2 (r)7 Q , / g(q, r) = h a p{(\, t{r)). For later convenience we fix the Lagrangian coordinates 
q a to have physical dimension of length and the conformal time variable r to have dimension 
of time. As a consequence the spatial metric 7 a/ 3 is dimensionless, as is the scale-factor a(r) 
which must be determined by solving the Friedmann equations for a perfect fluid of dust 

a'\ 2 8nG 22 

Q b a — kc , (17) 

2— - ( — ) + kc 2 = . (18) 
a \a ) 

Here primes denote differentiation with respect to the conformal time r and k represents the 
curvature parameter of FRW models, which, because of our choice of dimensions, cannot be 
normalized as usual. So, for an Einstein-de Sitter universe n = 0, but for a closed (open) 
model we simply have k > (k < 0). Let us also note that the curvature parameter is related 
to a Newtonian squared time-scale kn through kn = kc 2 (e.g. Peebles 1980; Coles & Lucchin 
1995); in other words k is an intrinsically PN quantity. 

By subtracting the isotropic Hubble-flow, we introduce a peculiar velocity-gradient tensor 

= a cu a ; , - -b^ = V 7 7 7 // , (19) 

where u a = (1/a, 0,0,0). 

Thanks to the introduction of this tensor we can rewrite the Einstein's equations in a more 
cosmologically convenient form. The energy constraint becomes 



a" 



d 2 - ^ v d\ + + c 2 (n - Qk) = lQirGa 2 g b 5 , (20) 
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where lZ a p (7) = a 2 ^R a p{h) is the conformal Ricci curvature of the three-space, i.e. that 
corresponding to the metric 7^; for the background FRW solution = (1 + f q 2 )~ 2 8 a/ 3, 

one has TZ a (5 (^ FRW ) = 2/t<5 a /3 . We also introduced the density contrast 5 = (g — g b )/ 1 g b . 
The momentum constraint reads 



nila = ^ , (21) 

where the double vertical bars denote covariant derivatives in the three-space with metric 7^. 

Finally, after replacing the density from the energy constraint and subtracting the back- 
ground contribution, the extrinsic curvature evolution equation becomes 



. (22) 



The Raychaudhuri equation for the evolution of the peculiar volume-expansion scalar 1? 
becomes 

■&' + -■& + + 4TvGa 2 g b 5 = . (23) 

The main advantage of this formalism is that there is only one dimensionless (tensor) variable 
in the equations, namely the spatial metric tensor 7^, which is present with its partial time 
derivatives through t?°o [Eq.(19) above], and with its spatial gradients through the spatial Ricci 
curvature T^ a ^- The only remaining variable is the density contrast which can be written in 
the form 

5(q,r) = (l + 5o(q))[7(q,r)/ 7 o(q)]" 1/2 -l , (24) 

where 7 = det 7^/3. A relevant advantage of having a single tensorial variable, for our purposes, 
is that there can be no extra powers of c hidden in the definition of different quantities. 



2.2 Fluid— flow approach 

Following the fluid-flow approach, described in the classical review by Ellis (1971) [see also 
Ehlers (1993)], we can alternatively describe our system in terms of fluid properties, in our 
case matter density, volume-expansion scalar and shear tensor, and two geometric tensors, the 
electric and magnetic parts of the Weyl tensor defined above. The derivation of the equations 
reported below is thoroughly described by Ellis (1971) and will not be reported here. 

For most cosmological purposes it is convenient to adopt the conformal rescaling and FRW 
background subtraction described in the previous sub-section. Therefore, we can start by 
writing the continuity equation directly in terms of the density contrast 5, 

^ + (1 + 5)^ = 0, (25) 

with -j^p denoting convective differentiation with respect to the conformal time r. In our 
Lagrangian frame, however, and for a scalar field, convective differentiation and partial dif- 
ferentiation coincide. The formal solution of this equation is given by Eq.(24) above. The 



10 



peculiar volume-expansion scalar d obeys the Raychaudhuri equation which we can rewrite in 
the form 

— + ^ + -$ 2 + a a ^ a + AiiGa 2 Q b 5 = , (26) 



Dt a 3 

P = U P ~ 3 U P^ 



where a a = "& a — \5 a 'd is the shear tensor. The shear, in turn, evolves according to 



Da a n' 2 1 

-pf + -°*p + 3^ + - 3*y v 7 + £a P = » ( 27 ) 

where we have rescaled the electric tide as = a 2 -E°g, which can be written in terms of our 
new variables as 

= ^> v ; + \^ a p + a -°% - + c 2 (n - l -5%n) . (28) 

Note that, for a generic second rank tensor A a p, one has 

Dr ~ dr 1 13 p 1 ' [ ' 

where 4- denotes the total derivative with respect to r, which in comoving coordinates coincides 
with the partial one. The two last terms in the r.h.s. come from writing the Christoffel 
symbols in our gauge. It is then clear that when the operator acts on either the shear 
or the complete i?°Jg tensor, the second and third term in the r.h.s. cancel each other and 
the convective and total differentiation coincide. This cancellation also occurs for a generic 
A a p if either the relevant Christoffel symbols vanish (as it is the case for the Newtonian limit 
in Eulerian coordinates) or the convective derivative acts on the eigenvalues of A a o and such 
a tensor has the same eigenvectors of a a p [as it is the case for the electric tide in the "silent 
universe" case (Barnes & Rowlingson 1989; Matarrese, Pantano & Saez 1993; Bruni, Matarrese 
& Pantano 1995b)]. 

The electric tidal tensor evolves according to 

D£ a n' 3 / 

P i _p<* i fife* , ra 7 fS __( ca 1 , a cl 

Dt + a ^ P + 13 + fit 7 2 y 7 P 7 P 



°1 

2 



7^(^ 7 ^ a 7 || 5 + T lS ^\\\)j + ^Ga 2 Qb {\ + 6)a a p = , (30) 



where we have rescaled the magnetic tide as 7i. a p = a 2 H a p and redefined the Levi-Civita tensor 
so that r] a/3 ' y = <y~ 1 / 2 e a/37 (for simplicity we used the same symbol after rescaling). 
Finally, the magnetic Weyl tensor evolves according to 



Dt 



+^Pv(^ s £% ll 8 + v a ' tS £\ ll s)=o. (31) 
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Note that, following the discussion above, in the last two equations the convective time deriva- 
tive must include the two terms proportional to the shear, as in Eq.(29). Note that, apart from 
the cases listed after Eq.(29), these two terms cannot be disregarded even in the Newtonian 
limit. 

In the fluid-flow approach, besides the evolution equations, one has to satisfy a number of 
constraint equations. One has: the momentum constraint, which we rewrite in the form 

= . (32) 

the Ji-a constraint (which we actually used in Section 2 to define the magnetic tide in terms 
of derivatives of the spatial metric), 

n% = (irvw + irv,,,,) , (33) 

the div 8 constraint, 

= -1^1^ o\H^ + ^-a 2 Qb 5,p , (34) 

and the div H constraint 

c 2 n a ma = lp M ^a\£^ . (35) 

In the above equations one also needs to know the three-metric 7 a/ g. This can be obtained 
from the evolution equation 

7a/j' = lloc^p , (36) 

which is however only valid in our Lagrangian coordinates. In order to completely fix the 
spatial dependence of the metric one also needs to specify the energy constraint (the trace of 
the Gauss-Codacci relations), which we rewrite in the form 

c 2 (U - 6k) = a a y a - -tf 2 - 4-# + 16nGa 2 Q b 5 . (37) 

3 a 

Although we will not use the fluid-flow approach in this paper it is interesting to have the 
complete form of the equations, with the correct powers of c 2 included, in order to understand 
the Newtonian meaning of the electric and magnetic tide. We will come back to this point in 
Section 4.1. 

2.3 Local Eulerian coordinates 

Our intuitive notion of Eulerian coordinates, involving a universal absolute time and globally 
flat spatial coordinates is intimately Newtonian; nevertheless it is possible to construct a local 
coordinates system which reproduces this picture for a suitable set of observers. This issue 
has been already addressed by Matarrese, Pantano & Saez (1994a,b), who introduced local 
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Eulerian - FRW comoving - coordinates x A which are related to the Lagrangian ones q a via 
the Jacobian matrix with elements 

J A M T) = — =8 A a + Pl(q, r), A = 1, 2, 3 , (38) 

where T> A a (q, r) is called deformation tensor. Each matrix element J A a labelled by the Eulerian 
index A can be thought as a three-vector, namely a triad, defined on the hypersurfaces of 
constant conformal time. As shown in (Matarrese, Pantano & Saez 1994a,b), they evolve 
according to 

J A d = ®\J\ , (39) 

which also follows from the condition of parallel transport of the triads relative to q along the 
world-line of the corresponding fluid element D(aJ' A a )/ Dt = (see also Kasai 1995). 

Our local Eulerian coordinates are such that the spatial metric takes the Euclidean form 
Sab, i-e. 

7a/3 (q, r) = 5 AB J A a (q, r)J B p(q, r) . (40) 
Correspondingly the matter density can be rewritten in the suggestive form 

f?(q,r) = ft(r)(l + 5o(q))[J(q,r)/J (qr 1 , (41) 

where J = detJ" A a . Note that, contrary to the Newtonian case, it is generally impossible in 
GR to fix Jq — 1, as this would imply that the initial Lagrangian space is conformally flat, 
which is only possible if the initial perturbations vanish. 



2.4 Linear perturbation theory in Lagrangian coordinates 

In this subsection we will deal with the linearization of the equations obtained in Section 2.1. 
This will be done mostly for pedagogical purposes, in that it will allow us to obtain a number 
of results which will turn out to be useful for the 1/c 2 expansion. Apart from this, it can be 
interesting to re-obtain the classical results of linear theory in the comoving and synchronous 
gauge only in terms of the spatial metric coefficients. 

Let us then write the spatial metric tensor of the physical (i.e. perturbed) space-time in 
the form 

la(S = la/3 + Waf3 , (42) 

with 7 Q/ 3 the spatial metric of the background space - in our case the maximally symmetric 
FRW one, *f a p = ia^ W ~ an d a small perturbation. Also, we assume that the only non- 
geometric quantity in our equations, namely the initial density contrast So, is everywhere much 
smaller than unity. 

As usual, we can take advantage of the maximal symmetry of the background FRW spatial 
sections to classify metric perturbations as scalars, vectors and tensors (e.g. Bardeen 1980). 
We then write 

W a p = Xlo.fi + C|a/3 + ~(£a|/3 + f/3|a) + , (43) 
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with 

e la = K a a = ^\ a = 0, (44) 

where a single vertical bar is used for covariant differentiation in the background three-space 
with metric 7^. In the above decomposition x an d C represent scalar modes, £ a vector modes 
and TT a tensor modes (indices being raised by the contravariant background three-metric). 

Before entering into the discussion of the equations for these perturbation modes, let us 
quote a result which will be also useful in the next sections. In the fi^ evolution equation and 
in the energy constraint the combination V a p 
first order in the metric perturbation one has 



4:71% — (1Z + 2k) 5% and its trace appear. To 



V%w) = -2 



(V 2 - 2/6)7^ + x\ a p + KxS a /> 



(45) 



where V 2 (-) = (-)| 7 7 - Only the scalar mode x an d the tensor modes contribute to the three- 
dimensional Ricci curvature. 

As well known, in linear theory scalar, vector and tensor modes are independent. The 
equation of motion for the tensor modes is obtained by linearizing the traceless part of the i? ^ 
evolution equation, Eq.(22). One has 



7r Q/3 " + 2-vrV " c2 (V 2 - 2/€)7T a/3 = , 

(X 



(46) 



which is the equation for the free propagation of gravitational waves in a FRW background 
(compare with Eq.(3) in the Einstein-de Sitter case). The general solution of this equation is 
well known (e.g. Weinberg 1972) and will not be reported here. 

At the linear level, in the irrotational case, the two vector modes represent gauge modes 
which can be set to zero, £ a = 0. 

The two scalar modes are linked together through the momentum constraint, which leads 
to the relation 

X = Xo + k(C " Co) • (47) 

The energy constraint gives 



(V 2 + 3k) 



-C + (4irGa 2 Q b -Kc 2 )((-( )-c 2 Xo 



= SnGa 2 Qb5o , 



(48) 



while the evolution equation gives 



C" + 2 V = c\ ■ 



(49) 



An evolution equation only for the scalar mode ( can be obtained by combining together 
the evolution equation and the energy constraint; it reads 



(V 2 + 3k) 



C" + -C'-^Ga 2 g b (C-Co) 
a 



-8nGa 2 g b 5 . 



(50) 
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On the other hand, linearizing the solution of the continuity equation, Eq.(24), gives 

5 = 5o-^(V 2 + 3 K )(C-Co), (51) 
which replaced in the previous equation gives 

5" + -5' - A7iGa 2 g b 5 = . (52) 
a 

This is the well-known equation for linear density fluctuations, whose general solution can be 
found in (Peebles 1980). Once 5(r) is known, one can easily obtain ( and x, which completely 
solves the linear problem. 

Eq.(50) above has been obtained in whole generality; we could have used instead the 
well-known residual gauge ambiguity of the synchronous coordinates to simplify its form. In 
fact, C is determined up to a space-dependent scalar, which would neither contribute to the 
spatial curvature, nor to the velocity-gradient tensor. For instance, we could fix ( so that 
(V 2 + 3/t)Co = — 2(5 , so that the ( evolution equation takes the same form as that for S. 

In order to better understand the physical meaning of the two scalar modes x an d C let us 
consider the simplest case of an Einstein-de Sitter background (k = 0), for which a(r) oc t 2 . 
By fixing the gauge so that V 2 Co = — 25o we obtain x( T ) — Xo and 

Or) = ^Xor 2 + Sor- 3 , (53) 

where the amplitude B of the decaying mode is an arbitrary function of the spatial coordinates. 
Consistency with the Newtonian limit suggests Xo = — ^^o, with tp the initial peculiar 
gravitational potential, related to 5 through V 2 y?o = 4:TTGalg 0b 5 - We can then write 

C(r) = -T^or 2 + B r- 3 . (54) 
This result clearly shows that, at the Newtonian level, the linearized metric is 

la(S = 8af5 + C\a/3 , (55) 

while the perturbation mode x is already PN. Note that also the tensor modes are at least 
PN. 

These results also confirm the above conclusion that in the general GR case the initial 
Lagrangian spatial metric cannot be flat, i.e. Jq ^ 1, because of the initial "seed" PN metric 
perturbation xo- 

3 Newtonian approximation 

The Newtonian equations in Lagrangian form can be obtained from the full GR equations of 
Section 2.1 by an expansion in inverse powers of the speed of light; as a consequence of our 
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gauge choice, however, no odd powers of c appear in the equations, which implies that the 
expansion parameter can be taken to be 1/c 2 . The physical meaning of this expansion has 
been already outlined in Section 1. 

Let us then expand the spatial metric in a form analogous to that used in our linear 
perturbation analysis of Section 2.4: 

7«tf = 7«tf + ^£^ + o(^), (56) 

where we made explicit the c dependence of the metric perturbation. The actual convergence 
of the series requires that the PN metric perturbation j^w^^ is much smaller than the 
background Newtonian metric 7 Q| g. Let us first concentrate on the Newtonian metric; the 
properties of w a p will be instead considered in Section 4. 

To lowest order in our expansion, the extrinsic curvature evolution equation, Eq.(22), and 
the energy constraint, Eq.(20), imply that V a p = V a p(i) = 0, and recalling that k = k, n /c 2 , 

K% = K a ^) = : (57) 

in the Newtonian limit the spatial curvature identically vanishes (e.g. Ellis 1971). This im- 
portant conclusion implies that j a/ 3 can be transformed to 5ab globally, i.e. that one can 
write 

lap = 5 AB J A a J% , (58) 

with integrable Jacobian matrix coefficients. In other words, at each time r there exist global 
Eulerian coordinates x A such that 

x(q,r) = q+S(q,r), (59) 



where S(q, r) is called the displacement vector, and the deformation tensor becomes in this 
limit 



dq a 

The Newtonian Lagrangian metric can therefore be written in the form 

/ \ s ( sA 9S A {q,r)\f^ B dS B {q,r)\ 
7c*(q,T) = 5 AB [5 A a + ™ [5% + ™ ' ) . (61) 

We can rephrase the above result as follows: the Lagrangian spatial metric in the Newtonian 
limit is that of Euclidean three-space in time-dependent curvilinear coordinates q a , defined at 
each time r in terms of the Eulerian ones x A by inverting Eq.(59) above. As a consequence, 
the Christoffel symbols involved in spatial covariant derivatives (which we will indicate by a 
single bar or by a nabla operator followed by greek indices) do not vanish, but the vanishing 
of the spatial curvature implies that these covariant derivatives always commute. 

Contrary to the evolution equation and the energy constraint, the Raychaudhuri equation, 
Eq.(23) and the momentum constraint, Eq.(21), contain no explicit powers of c, and therefore 
preserve their form in going to the Newtonian limit. These equations therefore determine the 
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background Newtonian metric j a p, i.e. they govern the evolution of the displacement vector 
S. 

The Raychaudhuri equation becomes the master equation for the Newtonian evolution; it 
takes the form 

& + + + AixGa 2 6b (T 1/2 - 1) = , (62) 

(X 

where 

^ = \l ai % P , (63) 

and, for simplicity, we assumed So = (a restriction which is, however, not at all mandatory). 
We also used the residual gauge freedom of our coordinate system to set 7^(10) = 5 a p, implying 
jfo = 1) i- e - to make Lagrangian and Eulerian coordinates coincide at the initial time. That 
this choice is indeed possible in the Newtonian limit can be understood from our previous 
linear analysis, where this is achieved by taking, e.g., ( = 0. 
The momentum constraint, 

K\» = t» , (64) 

is actually related to the irrotationality assumption. We will come back to this point in the 
next section. 

Before closing this section, let us notice a general property of our expression for the La- 
grangian metric: at each time r it can be diagonalized by going to the local and instantaneous 
principal axes of the deformation tensor. Calling 7^ the eigenvalues of the metric tensor, J a 
those of the Jacobian and d a those of the deformation tensor, one has 

7«(q, r) = J"j(q, r) = [1 + d a (q, r)f . (65) 

In Section 3.2 below, the diagonal form of the metric tensor will be reconsidered in the 
frame of the Zel'dovich approximation. Beyond the mildly non-linear regime, where this 
approximation is consistently applied, diagonalizing the metric is in general, i.e. apart from 
specific initial configurations, of smaller practical use, because metric (and deformation) tensor, 
shear and tide generally have different eigenvectors. 

From this expression it becomes evident that, at shell-crossing, where some of the Jacobian 
eigenvalues go to zero, the related covariant metric eigenvalues just vanish. On the other 
hand, other quantities, like the matter density, the peculiar volume expansion scalar and some 
eigenvalues of the shear and tidal tensor will generally diverge at the location of the caustics 
(see Bruni, Matarrese & Pantano 1995b, for a discussion). This diverging behaviour makes 
the description of the system extremely involved after this event. Although dealing with this 
problem is far outside the aim of the present paper, let us just mention that a number of 
ways out are available. One can convolve the various dynamical variables by a suitable low- 
pass filter, either at the initial time, in order to postpone the occurrence of shell-crossing 
singularities (e.g. Coles, Melott & Shandarin 1993; Kofman et al. 1994), or at the time 
when they form, in order to smooth out the singular behaviour (e.g. Nusser & Dekel 1992, 
and references therein); alternatively one can abandon the perfect fluid picture and resort to 
a discrete point-like particle set, which automatically eliminates the possible occurrence of 



17 



caustics, at least for generic initial data. At this level, anyway, we prefer to take a conservative 
point of view and assume that the actual range of validity of our formalism is up to shell- 
crossing. 



3.1 Jacobian approach 

A more direct way to deal with the Lagrangian Newtonian equations is to write them in terms 
of the Jacobian matrix J A a . This approach is obviously related to the more usual ones in 
terms of the displacement vector S or in terms of the deformation tensor T> A a (Buchert 1989; 
Moutarde et al. 1991; Bouchet et al. 1992; Buchert 1992; Catelan 1995). The evolution 
equation has been explicitly written directly in terms of the Jacobian matrix by Buchert & 
Gotz (1987), Lachieze-Rey (1993) and Catelan (1995). 

In order to rewrite the Raychaudhuri equation in terms of the Jacobian matrix, we notice 
that 

*y = " - , (66) 

where we have introduced the inverse Jacobian matrix 

J", - & • (") 

where Eulerian indices are raised and lowered by the Kronecker symbol. To make explicit our 
notation, we just stress that elements of will be characterized by a greek (i.e. Lagrangian) 

index subscript, while elements of the inverse matrix J^- will be characterized by a greek index 
superscript. 

Replacing the latter identity into the Newtonian expression Eq.(62) yields 

JaJ\" + -J' 1 J' = ^Ga 2 Q h (l - J' 1 ) . (68) 

(X 

Note that this expression is, apart from the use of a different time variable, identical to 
Eq.(60), in Catelan (1995) [see also Appendix A in (Buchert 1989), and (Buchert 1992)]. 
We also notice that the parallel transport condition, Eq.(39), can be rewritten in the form 

n = t a j a ; . (69) 

This equation, together with i?°Jg = |7 a7 7 7 // gives the general relation 

JaJJ A (] = JAaJ A p ■ (70) 
Replacing these relations in the momentum constraint we obtain in whole generality 

J a AJ%\ la + J a A\\ a J%' = (J'' J'),, ■ (71) 

On the other hand, in the Newtonian limit we have 
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as it follows from the fact that S A a p = S A ^ a . Using this commutation property it is easy to 
verify that 

r^ 7 = J\J p n ■ (73) 

Thanks to the latter relation and to the well-known matrix identity TrlnJ = IndetJ, it is 
straightforward to verify that the momentum constraint in the Newtonian limit becomes an 
identity. It is then clear that Eq.(70) is more fundamental than the momentum constraint: 
it plays the role of an irrotationality condition written in Lagrangian space. This is of course 
equivalent to the standard form [compare with Eq.(59) in (Catelan 1995)] 

e a ^J A p J' Al = . (74) 

This equation, together with the Raychaudhuri equation above, Eq.(68), completely determines 
the Newtonian problem, in terms of either the Jacobian matrix, the deformation tensor or the 
displacement vector. 

The very fact that we have been able to recover the standard equations for the Newtonian 
approximation in the Lagrangian picture, by starting from the Lagrangian GR treatment and 
expanding in powers of 1/c 2 , should be considered as a further confirmation of the validity of 
our method. 



3.2 Zel'dovich approximation 

Having shown the equivalence of our method, in the Newtonian limit, with the standard one, 
it is now trivial to recover the Zel'dovich approximation (Zel'dovich 1970). This is obtained 
by expanding Eq.(68) and Eq.(70) to first order in the displacement vector. The result is 

x(q,r)=q + J D(r)V<l>o(q), (75) 

where only the growing mode solution D(t) of Eq.(52) has been considered, and we introduced 
the potential 3>o(q), such that V^$o — So/Do, where V 2 q is the standard (i.e. Euclidean) 
Laplacian in Lagrangian coordinates; more in general, at this perturbative order covariant and 
partial derivatives with respect to the q a coincide. The potential $o is easily related to the 
initial peculiar gravitational potential defined in Section 1, $o — — {^Go^Qq^D^^lpq. 

More interesting is to derive from the above expression the corresponding Zel'dovich metric. 
It reads 

7 f| L (q,r) = v(^ + J D(r)$ ,^( q ))(5 5 /3 + J D(r)<I. , 5 /3 ( q )) . (76) 

One can of course diagonalize this expression by going to the principal axes of the defor- 
mation tensor. Calling \ a the eigenvalues of the matrix $o, a p, one finds 

7f £L (q,r) = [l + /J(r)A Q (q)] 2 . (77) 

Note that, contrary to what has been commonly done so far in the literature, the metric 
tensor must be evaluated at second order in the displacement vector, in order to obtain back 
the correct Zel'dovich expressions for the dynamical variables (density, shear, etc ...). 
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The above diagonal form of the metric allows a straightforward calculation of all the relevant 
quantities. The well-known expression for the mass density is consistently recovered, 



d** 11 = QbUO- + DX*)- 1 . (78) 

a 

The peculiar velocity-gradient tensor has the same eigenframe of the metric; its eigenvalues 
read 

= ilk ' (79) 

By summing over a the latter expression we can obtain the peculiar volume-expansion scalar 
and then the shear eigenvalues 

zel _ D \ a 1^ D X a 



1 + D\ a 3^1 + DX a 



3?rri' (81) 



The electric tide comes out just proportional to the shear. Its eigenvalues read 

£ z a EL = -AirGa^ a z a EL . (82) 

These expressions for the shear and the tide completely agree with those obtained by 
Kofman & Pogosyan (1995) and Hui & Bertschinger (1995). The fact that metric, shear and 
tide have simultaneous eigenvectors shows that fluid elements in the Zel'dovich approximation 
actually evolve as in a "silent universe" (Matarrese, Pantano & Saez 1994a; Bruni Matarrese 
& Pantano 1995b), with no influence from the environment, except for that implicit in the 
self-consistency of the initial conditions. 

So far the Zel'dovich approximation has been obtained by first taking the Newtonian limit 
(c — >• oo) of the GR equations and then linearizing them with respect to the Newtonian 
displacement vector. One could also drop the first step and linearize the GR equations of 
Section 2.1 with respect to the local deformation tensor as introduced in Section 2.3; in such 
a case one would get a fully relativistic version of the Zel'dovich approximation. 

The latter problem has been already discussed a number of times by various authors. Un- 
fortunately, there has been a lot of misunderstanding on what the "relativistic Zel'dovich ap- 
proximation" should actually be. Most authors just deal with the GR version of the Zel'dovich 
solution, i.e. with the non-linear evolution of planar perturbations, which is a sub-case of the 
well-known exact solutions obtained by Szekeres (1975). Such an approach, however, does 
not allow to deal with the approximate non-linear behaviour of generic perturbations in a 
relativistic framework. 
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3.3 Lagrangian Bernoulli equation 

As we have demonstrated above, it is always possible, in the frame of the Newtonian approx- 
imation, to define a global Eulerian picture. This will be the picture of the fluid evolution 
as given by an observer that, at the point x = q + S(q, r) and at the time r observes the 
fluid moving with physical peculiar three-velocity v = dS/dr. From the point of view of a 
Lagrangian observer, who is comoving with the fluid, the Eulerian observer, which is located 
at constant x, is moving with three-velocity c?q(x, r)/dr = —v. 

The line-element characterizing the Newtonian approximation in the Eulerian frame is well 
known (e.g. Peebles 1980) 



ds 2 = a 2 (r) 



- ( 1 + 2<Pg ^ ,Tj ) Jdr 2 + 5 AB dx A dx B 



(83) 



with ip g the peculiar gravitational potential, determined by the mass distribution through the 
Eulerian Poisson equation, 

V^,(x, r) = AnGa 2 (r)g b (r)5^, r) , (84) 

where the Laplacian V^, as well as the nabla operator V, have their standard Euclidean 
meaning. The perturbation in the time-time component of the metric tensor here comes from 
the different proper time of the Eulerian and Lagrangian observers. As already noticed in 
the Introduction, this Newtonian line-element is different to that of Eq.(2), describing the so- 
called weak-field limit; the extra term —2(p g /c 2 , multiplying the spatial line-element, would 
in fact give rise to PN corrections in the equations of motion for the matter. 

It is now crucial to realize that all the dynamical equations obtained so far, being entirely 
expressed in terms of three-tensors, keep their form in going to the Eulerian picture, only 
provided the convective time derivatives of tensors of any rank (scalars, vectors and tensors) 
are modified as follows: 

D 9 dS 

_^_ + v . V , VS -. (85) 

This follows from the fact that, for the metric above, T° AB = T^ B = T BC = 0, which also 
obviously implies that covariant derivatives with respect to x A reduce to partial ones. 

The irrotationality assumption now has the obvious consequence that we can define an 
Eulerian velocity potential through 

v(x,r) = V$,(x,r). (86) 

The Newtonian peculiar velocity-gradient tensor then becomes 

= • (87) 

because of which the momentum constraint gets trivially satisfied and the magnetic Weyl 
tensor becomes identically zero in the Newtonian limit. 
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We can now write the Raychaudhuri equation for the Eulerian peculiar volume-expansion 
scalar and use the Poisson equation to get, as a first spatial integral, the Euler equation 

v ' + v • Vv + -v = -Vy? 9 . (88) 

(X 

This can be further integrated to give the Bernoulli equation 

< + -^ + ^(V$,) 2 = -^. (89) 

CI Zi 

On the other hand, by taking gradients of the Euler equation we can obtain an Eulerian 
evolution equation for the tensor $ab- More interesting is that this equation can be transported 
back to the Lagrangian frame to get 

*V + + = -vf V (90) 

where must be thought as a Lagrangian peculiar gravitational potential to be determined 
through the Lagrangian Poisson equation 

<pf\ = 4nGa*g b (T 1/2 ~ 1) • (91) 

These two Lagrangian expressions will turn out to be very useful for the PN calculations 
of Section 4. 

There is, however, another consequence of these equations that we can easily derive. While 
in the Lagrangian frame the three- velocity field does not exist, the tensor t? q ^ is well-defined, 
so that we can rewrite Eq.(87) in Lagrangian coordinates to obtain a Lagrangian velocity 
potential <&( L \q, r), through 

^ = *§ )a p ■ (92) 

This potential obeys what we can name the Lagrangian Bernoulli equation, which is easily ob- 
tained from the Bernoulli equation above, provided we recollect the convective time derivative 
and express it in Lagrangian form. We get 

^ y + -$[ L) - V*<2*S = M L) • (93) 

The most astonishing difference between the Eulerian and Lagrangian versions of the Bernoulli 
equation is the relative sign of the temporal and spatial derivatives. We could obtain more 
similar forms by reversing the arrow of time and the sign of the gravitational interaction. In 
this sense, therefore, the Lagrangian Bernoulli equation acts as a sort of time machine (cf. 
Nusser & Dekel 1992). This fact becomes more clear if we think to the fact that, by solving 
it, we are indeed asking how the Lagrangian (i.e. initial) geometry at q should modify itself 
in order to reproduce the Eulerian (i.e. evolved) properties of the velocity and density fields 
at the point x(q, r) as time goes on. 
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This equation could be used in principle as an alternative Lagrangian formulation of Newto- 
nian theory, whose fundamental variables would be the velocity potential $^ L ) , the gravitational 
potential and the metric tensor j a p. This approach could be useful, in particular, in order 
to obtain new self-consistent approximation schemes to the non-linear evolution of dust in 
the Lagrangian frame. To this aim, however, we need two more equations to close the system. 
These can be provided by the Lagrangian Poisson equation above, Eq.(91), and by the very 
definition of d a p = t^t? 7 ^ which implies 

<C = \l' aP • (94) 

Of course, the Lagrangian scalars $^ and are related to their Eulerian counterparts 
by a simple coordinate transformation, namely <£>( L )(q, r) = $[, E - ) (x(q, r), r) and (p^ L \q, r) = 
^f)(x(q,r),r). 



4 Post— Newtonian approximation 

Having examined all the aspects of our formalism in the Newtonian limit, we are now ready to 
proceed to the next perturbative order in 1/c 2 . The PN terms ^w^' in Eq.(56) should be 
thought as small perturbations superposed on a Newtonian background *j a p. The fact that the 
three-metric in the Newtonian limit is that of Euclidean space in time-dependent curvilinear 
coordinates q a , implies that we can apply most of the standard tools of linear perturbation 
theory in a flat spatial background. In particular, we can classify our PN metric perturbations 
as scalar, vector and tensor modes, as usual. 
We then write 

,„(™0 - „( PN )~ 4- r (PN) 4- 1 (F {PN) 4 ^ (PN h 4 J PN) (Q^\ 
W af3 -X" 7a/3 + C| a/3 + T^l/ 3 + S 9 !" > + a P ' ^ > 

with 

Z (PN)a la = - (PN)a a = - {PN)a , la = 0, (96) 

where greek indices after a single vertical bar, or nabla operators with a greek index, denote 
covariant differentiation in the Newtonian background three-space with metric 7 aj g. In the 
above decomposition x^ PN ^ an d represent PN scalar modes, PN vector modes and 

tt^ N ^ PN tensor ones (indices being raised by the contravariant background three-metric). 
We deliberately used the same symbols as in Section 2.4, in order to emphasize the analogy 
with the linear problem. Some of these PN modes, namely x^ PN ^ an d ^% N \ also have a non- 
vanishing linear counterpart, as noticed in Section 2.4 (actually the linear part of appears 
as a gauge mode in the equations), while others, namely C,^ PN \ and £^ PAr) are intrinsically 
non-linear. Unlike linear perturbation theory in a FRW background, metric perturbations of 
different rank do not decouple: this is because our time-dependent Newtonian background 
enters the equations not only through the metric 7 Q/ g, but also through the peculiar velocity- 
gradient tensor & a 0: which also contains scalar, vector and tensor modes. This fact leads to 
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non-linear scalar-vector, scalar-tensor and vector-tensor mode mixing, which also explains 
why we had to account for the vector modes ^ PJV - ) in the expansion of w^ N \ in spite of the 
irrotational character of our fluid motions. Actually, that vector modes appear in the non- 
linear evolution of an irrotational fluid in Lagrangian coordinates is well known also in the 
Newtonian framework (see Buchert 1994; Catelan 1995). 

As in every perturbative calculation, some of the equations have the property to mix 
different perturbative orders. This is of course necessary in order to make the n-th order 
coefficients of the expansion calculable in terms of those of order n — 1. In our case the 
energy constraint and the extrinsic curvature evolution equation (which at the Newtonian 
level implies lZ a l3 ( , y) = 0) play this role. Therefore we assume that the Newtonian metric and 
its derivatives are known by solving the Raychaudhuri equation and the momentum constraint, 
and we calculate the PN metric perturbations in terms of them. 

Let us first compute the tensor V a p = 47?°^ — (1Z + 2k,)5°' /3 to first order in 1/c 2 . We obtain 

c*V%w^) = -2 ( VV™) Q , + , (97) 

where now V 2 (-) = (•){*■ 

In the energy constraint only the scalar \ enters, 

V 2 x = \ (> ~ KK) + 2^ - ^Ga 2 g b ( T 1/2 ~ 1) , (98) 

where, here and from now on, we have dropped the superscript (PN) on PN terms. One 
can also obtain an equation for \ from the trace of the evolution equation, Eq.(22), which is 
however equivalent to the latter, thanks to the Newtonian Raychaudhuri equation. 

The tensor perturbations are instead determined via the evolution equation, Eq.(22) 
(actually from its trace-free part), 

V 2 ^ = 20°; + 2(2^ + - V^ 2 - V,) -x\%- (99) 

A by-product of the latter equation is that linear tensor modes, which in the c — > 00 limit 
appear as harmonic functions (i.e. pure gauge modes), do not contribute to the r.h.s., i.e. to 
the Newtonian evolution of the system, as expected. 

In order to get an equation for the tensor modes decoupled from the scalar mode x we 
can resort to the equations obtained in Section 3.3 above. To this aim we define the auxiliary 
Newtonian potential ^[ L \ through 

vM L > = ~ - v M ) • (100) 

Using this definition in Eq.(98), we obtain 

X = -*i L) + 2-<t>l L) - 2^> (101) 

a a 
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(which, at the linear level reduces to the expression of Section 2.4, xo = _ ^Vo)- Using this 
expression in Eq.(99) and replacing i? ^' from Eq.(90) we obtain 

V 2 tt^ = (V V/3 + 5^ V 2 ) + 2 - tf^) , (102) 

which has the significant advantage of being explicitly second order (in any possible perturba- 
tive approach). This equation is one of the most important results of this paper: it gives (in 
the so-called near zone) the amount of gravitational waves emitted by non-linear cosmological 
perturbations, evolved within Newtonian gravity. In other terms, this equation, which is only 
applicable on scales well inside the horizon, describes gravitational waves produced by an inho- 
mogeneous Newtonian background. At first sight it may appear surprising that gravitational 
waves already appear in a PN calculation, whilst in the standard calculations in non-comoving 
gauges one finds them at the PPN level. This is indeed a peculiarity (actually an advantage) 
of the Lagrangian coordinates, where all orders are scaled by two powers of c: for instance, 
the Newtonian terms already appear at the zeroth order, whereas in the longitudinal gauge 
the gravitational potential carries a 1/c 2 factor. 

This formula can be compared with the PN limit of Eq.(3) in Section 1, to which it actually 
reduces (up to an already mentioned numerical factor) if the i5 a o are calculated from linear 
theory and in an Einstein-de Sitter model. We have ^"^(q, r) = D'(r)<&o"g(q), with D(r) the 
growing mode solution of Eq.(52) (D(t) oc a(r) oc r 2 in the Einstein-de Sitter case) and $o 
defined in Section 3.2. Therefore 



Y72 _a n /2 



^ + ^V^ + 2 $ ^V 2 $ - %, a ^o,^ 



(103) 



where the symbol V 2 indicates the standard (Euclidean) form of the Laplacian in Lagrangian 
coordinates and \l/ = ^v(tq) and indices are raised by the Kronecker symbol. 

To completely determine the PN metric perturbations we still need the scalar mode ( 
and the vector modes £ Q , which can be computed through the momentum constraint and the 
Raychaudhuri equation. We then need $ a 3 to PN order; it reads 

n = n + \ 2 ($ a y P - ^> a 7 + \<) ■ ( l04 ) 

By replacing this expression into the momentum constraint, Eq.(21), we obtain 

— A is- _ _ rT , 

-v,P 



-4«jv^ l 2, (105) 



where the term on the r.h.s. comes from the relation i?% a — $,p = 2/c<&^, which can be easily 

(L) a (L) a (T\ 9 

derived by expanding the Jacobi identity ^ — ^ = ^^TZ a s , in powers of 1/c . 

In order to write the last equation in terms of the various PN perturbation modes, the 
Newtonian identity T^ p = is also useful, implying 

K/)|a = K/3|a)' - KvTfi + *V>° ■ ( 106 ) 



25 



By replacing the expansion of w a p into this equation we finally get 



-^| a + - *W« + ^V 2 e« = 4k n ^\ , (107) 

which, at the linear level, reduces to Eq.(47) above. 

One can verify that the three-divergence of the last equation reduces to an indentity, 
therefore, in order to completely determine the three remaining PN modes ( and £ a we need 
one more equation. This is in fact provided by the PN Raychaudhuri equation, which reads 

w" + -w' + 2#>V = 4nGa 2 g b {l + 5){w- w ) , (108) 

(X 

having assumed Sq PN ^ = 0, as suggested by linear theory. By replacing the expansion of w a p 
into this equation we get 

(3 X + V 2 C)" + ^(3 X + V 2 C)' + 2$ X ' + + = 



= AnGa 2 Q b {l + 5) 



3(x - Xo) + V 2 C 



(109) 



where we have also set w = 3xo, in agreement with linear theory. 

Unfortunately, we have not been able to further simplify these equations, which nevertheless 
show that ( and £ Q are implicitly determined by the Newtonian quantities, once \ an d ^% 
have been computed. There is a caveat concerning second-order perturbation theory, where 
( is left undetermined by the above equations. Nevertheless, such an ambiguity is completely 
removed by going to the PPN energy constraint, which can be easily shown to fix w and then 
( without explicit knowledge of the truly PPN coefficients. Let us complete this discussion 
by reporting the PPN energy constraint. Defining a PPN metric perturbation w^ PPN ^ /c 4 , one 
gets 

(tf + 2 -)w' - + = -8nGa 2 g b (l + S)(w - w ) , (110) 

with 7?/ PPAr ) the PPN conformal Ricci scalar, 

tf pp * = -2k n w + w( PPN ^ Vv - V V PPJV ) + vfM\ + - 2 ^ 7 \J + 

+l<i r ^ - l w V w ^ - - 2 <i,)(^r - 2^, 7 ) . (in) 
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4.1 Fluid— flow approach in the Newtonian limit 

We are now ready to discuss the fluid-flow approach presented in Section 2.2, within the 
Newtonian approximation. The reason why this discussion has been included in this Section is 
that, as we shall see, some of the relevant tensors must be computed at the PN order in order 
to provide the correct Newtonian evolution of the system. 

We just have to discuss the order in our 1/c 2 expansion at which the various tensors 
enter the equations of Section 2.2. It is immediately clear that the mass continuity equation, 
Eq.(25), the Raychaudhuri equation, Eq.(26), and the shear evolution equation, Eq.(27), where 
no explicit powers of c appear, just keep their form, once the various tensors are replaced by 
their Newtonian counterparts. So, we have 

S'+ (1 + 5)0 = 0, (112) 
& + + \$ 2 + + ^Ga 2 g b 5 = , (113) 

and 



< + + + a^o\ - -5^a\a\ + £ a p = . (114) 



On the other hand, by its very definition, Eq.(28), the electric tide contains a contribution 
coming from the PN terms x an d Tt a p because of the spatial curvature terms, 



1 _ U „. a' „. „. „ 1 



i " 3" u" p -r p -r —° (3 ~ & 7 ^ 7 /3 ~ 2 



V^+(V a V r ^V 2 ]x 



(115) 



It is however immediate to realize that, once the expressions of this section for the PN tensors 
X and ix a p are used, one recovers the simpler form, 

£*(,= (v a V^-^V 2 )yf), (116) 

which, in Eulerian coordinates reduces to the standard form 

£ab = ¥ 9 ,AB - ^ AB V 2 x ip g . (117) 

On the other hand, if we replace in Eq.(33), the Newtonian peculiar velocity-gradient ten- 
sor, we obtain the well-known result (e.g. Ellis 1971) that the magnetic tensor identically 
vanishes in the Newtonian limit. This can be very easily shown by either applying the formal- 
ism of Section 3.3, i.e. writing $ a g through covariant derivatives of the Lagrangian velocity 
potential, or by writing the same tensor in terms of the Jacobian matrix of Section 3.1. The 
physics underlying this result is the conformal flatness of the Newtonian spatial sections, im- 
plying the commutation of spatial covariant derivatives. A simple consequence of this fact is 
that, at the Newtonian level, the div £ constraint, Eq.(34), reduces to 
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which, owing to our expression for £, turns out to be just the gradient of the Lagrangian 
Poisson equation, Eq.(91), namely 

VVf ) = ^Ga 2 g h 5 . (119) 

Let us now come to the tide evolution equation, Eq.(30). In that evolution equation the 
circulation of the magnetic tensor is multiplied by c 2 , which means that the PN part of curl 7i 
is the source of non-locality in the Newtonian electric tide evolution equation. On the other 
hand, if we look at the magnetic tide evolution equation, which starts to be non-trivial at the 
PN order, we see that curl £ is consistently a PN quantity. 

The Newtonian meaning of the momentum constraint, Eq.(32), has been already discussed 
in Section 3. Also interesting is the div 7i constraint, Eq.(35), telling us that the general 
non-vanishing of the PN magnetic tensor (see also Lesame, Ellis & Dunsby 1996), implies that 
the Newtonian shear and electric tide do not commute, i.e. they have different eigenvectors 
(viceversa, their non-alinement causes a non-zero div 7i). Another possible version of this 
result is that the ratio of the velocity potential to the gravitational potential beyond the linear 
regime becomes space-dependent. 

To summarize our results, we can say that within the Newtonian approximation the fluid- 
flow approach in Lagrangian coordinates can be formulated in terms of mass continuity, Ray- 
chaudhuri and shear evolution equations plus the Newtonian div £ constraint, which closes 
the system, provided we remind the circulation-free character of the electric tide in this limit. 
Of course the direct use of a constraint to close the system of evolution equations, has the 
disadvantage of breaking the intrinsic hyperbolicity of the GR set of evolution equations, so 
that the entire method looses its basic feature. No ways out: this is the price to pay to the 
intrinsic non-causality of the Newtonian theory [see also Ellis (1990)]. 

The above discussion on the role of the PN magnetic tidal tensor, as causing non-locality 
in the Newtonian fluid-flow evolution equations, completely agrees with a similar analysis 
by Kofman & Pogosyan (1995). The only variant is that we obtained our results directly in 
Lagrangian space, while they worked in non-comoving (i.e. Eulerian) coordinates. A different 
point of view on the subject is expressed by Bertschinger & Hamilton (1994), according to 
which the magnetic part of the Weyl tensor is non-vanishing already at the Newtonian level. 
According to Kofman & Pogosyan (1995) the difference might be "semantic"; most important, 
there is general agreement on the fundamental fact that the Newtonian tide evolution is affected 
by non-local terms. 

4.2 Post— Newtonian tensor modes within a collapsing homogeneous 
ellipsoid 

The PN expression for TT a/ 3, Eq.(102), has the relevant feature of being non-local, through 
the presence of the scalar \EV A simpler way to deal with this problem is to transform the 
equation in Eulerian form, where it is easier to deal with the Laplacian operator V 2 (which 
has there the standard Euclidean form), obtain the Eulerian gravitational- wave tensor ixab 
and then go back to the Lagrangian expression through vr Q/3 (q, r) = J A a J b /3 ttab(^-((1, t), t). 
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We immediately obtain the Eulerian expression 

VI-kab = + 5 AB V 2 X ¥^ + 2 (Mm - # AC # C B ) , (120) 

with 

v^f) = -i(^-^n), (i2i) 

which generally allows a simpler derivation of ir AB , given the (gradients of the) velocity po- 
tential. For a general homogeneous and isotropic random field $„, for instance, tt ab can be 
obtained by a simple convolution in Fourier space. Nevertheless, we would like to obtain 
here an analytic estimate of this tensor, in some simple cases. What we need is a model for 
non-spherical and non-planar collapse. The simplest model we can figure out is that of a 
homogeneous ellipsoid with uniform internal overdensity S(t) with respect to a FRW back- 
ground, of density Qbij) and scale factor a(r), in which it is embedded (White & Silk 1979; 
Peebles 1980). Calling R a (t) = (i(t)X a (t), A = 1,2,3, the physical length of the three axes, 
the peculiar gravitational potential within the ellipsoid is given by the simple expression 

y? 9 (x, r) = nGa 2 g b S ^ a A x 2 A , (122) 

A 

The dimensionless structure constants a A are defined in terms of the three X A through 

roc 

a A = X,X 2 X 3 / ds{X\ + s)- 1 H(X 2 B + s)' 1 ' 2 (123) 

B 

and are normalized so that Y,A a A = 2. In the particular case of oblate or prolate spheroids 
one can get explicit expressions for the a A in terms of an eccentricity parameter (e.g. Kellogg 
1953). The simplest non-trivial case, however, is that of an infinite cylinder, for which the two 
non-vanishing structure constants have the value 1 at any time. The intrinsic self-similarity of 
the equations of motion for fluid elements within the object implies that the overall shape and 
homogeneity are preserved at all times. One then usually makes the reasonable approximation 
that also the universe outside the ellipsoid stays uniform. In such a case, the ellipsoid axes 
obey the equation 

a' 

X" A + —X' A = -2ixGa 2 Q b 5 a A X A , (124) 

(X 

while mass conservation implies 

(1 + 5)X 1 X 2 X 3 = const . (125) 

The peculiar velocity-gradient tensor has eigenvalues r d A = X' A jX A (which remain unchanged 
in Lagrangian coordinates). These determine the Eulerian potential through the equation 

VM E) = -EWa-i, (126) 

A 
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where we adopt a notation such that A — 1 = 3 if A = 1 and A+l — 1 if ^4 = 3. Accounting 
for the ellipsoidal symmetry of the problem, this equation is solved by 

*f ) = ~E^-iE«fl4- (127) 

A B 

Replacing this solution in Eq.(120) gives 

V 2 X 1TA = ^ A {t) = --(0A0A-1 + $A$A+1 + + ($A$A-1 + $A$A+1 ~ $A-l$A+l) , 

(128) 

where ha = t^aa (no summation over repeated indices is understood) indicates a diagonal 
component, and J2a^a = 0. This equation is solved by 

TTA = jfJ-A E a BX 2 B . (129) 
4 B 

The off-diagonal components are instead harmonic functions, 

VIttab = , A ^ B . (130) 

These equations must be solved accounting for the transversality condition J2a n AB,A = 0, 
which gives 

kab = -^abx a x b , A^ B (131) 
(no summation over repeated indices), with vab = ^ba and 



^12 


= At3«3 


- UxOti 


~ /i 2 «2 




= /i2«2 


- fliOli 


- /^3«3 


V<Fo ■■ 


= - 


- HlOLl - 


- yU3«3 



(132) 



These formulae are completely general and do not contain approximations, (apart from those 
implicit in the homogeneous ellipsoid model). The results of a numerical integration of 
Eqs.(123), (124) and (125) are shown in Figure 1, for the collapse of a triaxial homogeneous 
ellipsoid. 

One might also use them to get the gravitational- wave emission outside the object (i.e. in 
the wave zone), by suitable matching with the interior solution. Far away from the body, one 
could, of course, also obtain the emitted gravitational radiation in terms of the quadrupole 
moments of our homogeneous ellipsoid (e.g. Landau & Lifshitz 1980, Sects. 41 and 105) 
and recover the usual 1/c 5 dependence of the radiated energy. In the wave zone, moreover, 
the transverse and traceless character of the ttab would allow to define the two polarization 
states of the graviton. These problems will not be considered here; it is nevertheless important 
to realize that our formalism is completely consistent with the qualitative expectations of 
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the quadrupole approximation, as it implies identically vanishing tensor modes in the case of 
spherical (cci = a 2 = a 3 = 2/3) and plane-parallel (cci = a 2 = 0, a 3 = 2) collapse. 

What we are interested in here is the behaviour of the PN tensor modes when the object 
is close to collapse. Of course, the set formed by Eq.(123), Eq.(124) and Eq.(125) could be 
integrated numerically to get the time evolution for the axes X A , the eigenvalues $a, and the 
structure constants a a (e.g. White & Silk 1979). To catch the qualitative behaviour close to 
collapse, however, we can safely apply the Zel'dovich approximation, which, for the evolution 
of the axes, yields 

X a (t) = X a (t ){1 + D(t)X a ) , (133) 

with A a = — jDo~ a A( T o)- These expressions should then be replaced into the definition of the a a 
to get them self-consistently. Nevertheless, according to White & Silk (1979) a rough estimate 
is obtained by simply neglecting the time-dependence of the a a- One can immediately derive 
the Jacobian eigenvalues J a = 1+D\ a) with A = a, and those of the peculiar velocity-gradient 
tensor $ a = D'\ a /(\ + D\ a ). 

These expressions can then be replaced into the previous equations to get the Lagrangian 
relations 

7i a = 7T aa (q, r) = -Ai a Jl a ^hl ' ( 134 ) 
4 p 

for the diagonal components, and 

7Ta/3(q, T) = ^VapJaJpqaqp , <* ^ (3 , (135) 

for the off-diagonal ones, with 



D 
-J 



/2 r 



1 



Aa-iA a+1 + \ a \ a -i + A a A a+1 + 3D\ a \ a _i\ a+ i )a a + 

xi (136) 
+ ( A a A Q _i + A a A a+ i — A a _iA a+ i + D\ a \ a -i\ a +i I 

and h> a /3 calculated from these fi a according to Eq.(132). For the most typical case of pancake 
collapse, where one Jacobian eigenvalue goes to zero first, these expressions also go to zero, 
like J. 

At this point we are able to compare the behaviour of these PN tensor modes to that of 
the Newtonian part of the metric, which is diagonal with eigenvalues 7 Q = J^. It is then 
clear that these PN modes vanish more slowly than the Newtonian part; their ratio diverges 
like i.e. like the mass density at collapse. Using a more refined approximation for the 

axes evolution, such as the one proposed by White & Silk (1979) or that recently proposed 
by Hui & Bertschinger (1995), would not change this qualitative result. Indeed, the numerical 
integration we have performed shows that such a divergence can be even stronger, as clearly 
displayed by the last panel in Figure 1. 

The homogeneous ellipsoid model we have worked out does not allow, unfortunately, to 
distinguish the global collapse from a shell-crossing singularity, but we may argue that this 
qualitative behaviour would generally apply even at shell-crossing. 
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5 Conclusions 

In this paper we have considered a Lagrangian approach to the evolution of an irrotational and 
collisionless fluid in general relativity. The use of a synchronous and comoving gauge allowed 
to reduce the fundamental variables of the system to the six metric tensor components of the 
spatial hypersurface orthogonal to the flow lines. Our method was based on a standard 1/c 
expansion of the Einstein and continuity equations which led to a new, purely Lagrangian, 
derivation of the Newtonian approximation. One of the most important results in this respect 
is that we obtained a simple and transparent expression for the Lagrangian metric; exploiting 
the vanishing of the spatial curvature in the Newtonian limit we were able to write it in terms 
of the displacement vector S(q, r) = x(q, r) — q, from the Lagrangian coordinate q to the 
Eulerian one x of each fluid element, namely 

ds 2 = a 2 (r) 

The spatial metric is that of Euclidean space in time-dependent curvilinear coordinates, con- 
sistently with the intuitive notion of Lagrangian picture in the Newtonian limit. Read this 
way, the complicated equations of Newtonian gravity in the Lagrangian picture become much 
easier: one just has to deal with the spatial metric tensor and its derivatives. The involved 
matrices appearing in the standard formulation are nothing else than the covariant and con- 
travariant metric tensor and the spatial Christoffel symbols, appearing in covariant derivatives. 
Moreover, the fact that the spatial Ricci curvature vanishes in this limit has the great practical 
advantage that spatial covariant derivatives commute. 

Next, we considered the post-Newtonian corrections to the metric and wrote equations for 
them. In particular, we were able to derive a simple and general equation for gravitational- 
wave emission from non-linear structures described through Newtonian gravity. The result is 
expressed in Lagrangian coordinates by Eq.(102), but it can also be given the Eulerian form of 
Eq.(120). These formulae allow to calculate the amplitude of the gravitational-wave modes in 
terms of the velocity potential which in turn can be deduced from observational data on 
radial peculiar velocities of galaxies, applying the POTENT technique (Bertschinger & Dekel 
1989). 

In the standard case, where the cosmological perturbations form a homogeneous and 
isotropic random field, we can obtain a heuristic perturbative estimate of their amplitude 
in terms of the rms density contrast and of the ratio of the typical perturbation scale A to the 
Hubble radius r H = cH~ l . One simply has 




(138) 



as it can be easily deduced from Eq.(103), specialized to an Einstein-de Sitter model. This 
effect gives rise to a stochastic background of gravitational waves which gets a non-negligible 
amplitude in the so-called extremely-low-frequency band (e.g. Thorne 1995), around 10~ 14 - 



-c 2 dr 2 + 5ab 



da a 



/3 + 



da f: 



(137) 
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10 15 Hz. We can roughly estimate that the present-day closure density of this gravitational- 
wave background is 



Note that this background is mostly produced here and now, so its energy is not affected 
by the usual a -4 dilution of gravitational radiation within the Hubble radius. In standard 
scenarios for the formation of structure in the universe, the typical density contrast on scales 1 
- 10 Mpc implies that VL gw is about 10 -5 - 10~ 6 . We might speculate that such a background 
would give rise to secondary CMB anisotropies on intermediate angular scales: a sort of tensor 
Rees-Sciama effect. This issue will be considered in more detail elsewhere. 

On much smaller scales, where the effect might be even more relevant, pressure gradients 
and viscosity cannot be disregarded anymore and the entire formalism needs to be largely 
modified. 

However, our PN formula also applies to isolated structures, where the density contrast 
can be much higher than the rms value, and, what is most important here, shear anisotropies 
play a fundamental role, as it happens in the formation of pancakes. A calculation of 7r Q/ g 
in the simple case of a homogeneous ellipsoid showed that the PN tensor modes become 
dominant, compared to the Newtonian contributions to the metric tensor, during the late 
stages of collapse, and possibly even in the case of a shell-crossing singularity. There are a 
number of important limitations of this result, the most important of which is the role that 
pressure would certainly play during the highly non-linear stages. A possible consequence 
could be that pressure gradients halt the growth of anisotropy before our relativistic effects 
come into play. It is nevertheless important to stress that our effect generally contradicts the 
standard paradigm, according to which the smallest scale for the applicability of the Newtonian 
approximation is set by the Schwarzschild radius of the object. Such a critical scale is indeed 
only relevant for nearly spherical collapse, whereas our effect becomes important precisely if 
the collapsing structure strongly deviates from sphericity. On the other hand, if we consider 
the dynamics of a collisionless fluid as a formal problem on itself, the fact that PN terms 
dominate over Newtonian ones implies that in such a regime the perturbative 1/c expansion 
breaks down and one should resort to a fully relativistic approach. 
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Figure 1. Evolution of the PN tensor modes during the collapse of a triaxial homogeneous 
ellipsoid embedded in an Einstein-de Sitter background; as well known the generic triaxial 
case exhibits pancake collapse. The top left panel shows the behaviour of the physical axes of 
the ellipsoid, Ra(j) = o,{r)X A {j) (A = 1,2,3), vs. the FRW scale-factor a(r); the lengths of 
the three axes are initially scaled as 1 : 1.25 : 1.5; the initial density contrast is 5q = 0.05; the 
evolution of the structure constant a a is also shown {top right panel) . The bottom left panel 
shows the evolution of the quantities /ia [defined in Eq.(128)], in terms of which the PN tensor 
modes iiab are obtained; note that two of the \ia tend to diverge near the pancake collapse of 
the ellipsoid; as shown in the bottom right panel, such a divergence is even stronger than that 
of the comoving mass density, 1 + S(t). 



37 



